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Abstract 

Multiple polylogarithms appear in analytic calculations of higher order corrections in quan- 
tum field theory. In this article we study the numerical evaluation of multiple polylogarithms. 
We provide algorithms, which allow the evaluation for arbitrary complex arguments and 
without any restriction on the weight. We have implemented these algorithms with arbitrary 
precision arithmetic in C++ within the GiNaC framework. 



1 Introduction 



In recent years multiple polylogarithms [1-3] found the way into physics in the field of pertur- 
bative calculations in quantum field theory. They satisfy several useful algebraic properties, in 
particular they satisfy two different Hopf algebras. Multiple polylogarithms are generalisations 
of harmonic polylogarithms [4,5], Nielsen polylogarithms and the classical polylogarithms [6,7] 
to multiple scales. Polylogarithms have a natural grading defined by their weight. From explicit 
higher order calculations it is emerging that one can express the results of Feynman integrals 
in terms of multiple polylogarithms. In dimensional regularisation we obtain in the finite terms 
polylogarithms up to weight 21 for an /-loop amplitude. Whereas we expect harmonic polyloga- 
rithms to be sufficient to express the results for 2 — > 2 scattering processes in massless quantum 
field theories, amplitudes with more external particles and/or massive particles involve addi- 
tional scales, which naturally leads to multiple polylogarithms. Since many recent results [8-27] 
of higher order calculations are expressed in terms of multiple polylogarithms, a numerical eval- 
uation routine for multiple polylogarithms is of immediate use for perturbative calculations in 
particle physics. Algorithms for the numerical evaluation of multiple polylogarithms are the 
subject of this article. For a few specific subclasses of multiple polylogarithms numerical evalu- 
ation methods can be found in the literature: The numerical evaluation of Nielsen polylogarithms 
has been known for a long time [28-30]. A special case of multiple polylogarithms are multi- 
ple zeta values, which are obtained from the polylogarithms if all scales are equal to one. For 
this special case, efficient algorithms for the numerical evaluation have been studied by Cran- 
dall [31] and Borwein et al. [2]. Furthermore, two recent papers provide numerical routines 
for harmonic polylogarithms and two-dimensional harmonic polylogarithms [32-34]. However 
these last mentioned routines are restricted to not more than two scales and weight not higher 
than 4. Here we provide methods for the larger class of multiple polylogarithms without any 
restrictions on the weight and the number of scales. We have implemented these algorithms with 
arbitrary precision arithmetic in C++ within the GiNaC [35] framework. 

This paper is organised as follows: In Section[2]we review as an introductory example the nu- 
merical evaluation of the dilogarithm. SectionEJdefines the notation for multiple polylogarithms 
used in this paper. Section|3]collects several useful properties of multiple polylogarithms, which 
will be used to construct the algorithms for the numerical evaluation. Section |5] is the main part 
of this paper and gives the algorithms for the numerical evaluation. In Sectional we report on 
the implementation in C++ within the GiNaC framework and describe the checks that we have 
performed. Finally Section |7] contains our conclusions. 

2 The dilogarithm 

In this short section we review as an introductory example the numerical evaluation of the dilog- 
arithm [30]. The numerical evaluation algorithm is rather simple, but does contain many ideas 
which will reappear in more elaborate form in the remaining part of the paper. The dilogarithm 
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is defined by 



.V 



, ln(l-f) 

U 2 (jc) = - / A— -, (1) 



and has a branch cut along the positive real axis, starting at the point x = 1. For |x| < 1 one has 
the convergent power series expansion 



00 x n 



Li 2 « = (2) 

n=\ 

The first step for a numerical evaluation consists in mapping an arbitrary argument x into the 
region, where the power series in eq. © converges. This can be done with the help of the 
reflection identity 

Li 2 (x) = -Li 2 Q)-^-i(ln(-*)) 2 , (3) 

which is used to map the argument x, lying outside the unit circle into the unit circle. The function 
ln(— x) appearing on the r.h.s. of eq. © is considered to be "simpler", e.g. it is assumed that a 
numerical evaluation routine for this function is known. In addition we can shift the argument 
into the range — 1 < Re(x) < 1/2 with the help of 

Li 2 (x) = -Li 2 (l-x) + — -ln(x)ln(l-x). (4) 

6 

Although one can now attempt a brute force evaluation of the power series in eq. ©, it is more 
efficient to rewrite the dilogarithm as a series involving the Bernoulli numbers B n : 

with z = — ln(l — x). The Bernoulli numbers B n are defined through the generating function 

+ °° fit 

~r~r = L B n~- (6) 

e — 1 ~ n n\ 

Therefore the numerical evaluation of the dilogarithm consists in using eqs. © and © to map 
any argument x into the unit circle with the additional condition Re(x) < 1/2. One then uses the 
series expansion in terms of Bernoulli numbers eq. ©. 

It occurs quite frequently in physics, that the numerical value for the dilogarithm is needed 
for the case, where the variable x is real except for a small imaginary part. The small imaginary 
part specifies if the dilogarithm should be evaluated above or below the cut in the case x > 1 . In 
this case the evaluation is split into the evaluation for the real part and the one for the imaginary 
part. The imaginary part of the dilogarithm is related to the logarithm: 

L12 (x) = Re L12 (x) — i In (x) Imln ( 1 — x) . (7) 
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The imaginary part of the logarithm is given for real x by 

ln(-x=F*0) = ln(|x|)=F«te(x), (8) 

where the step function 0(x) is defined as 0(x) = 1 for x > and 0(x) =0 otherwise and "iO" 
denotes a small imaginary part. 

3 Definitions 

3.1 Definition of multiple polylogarithms 

Multiple polylogarithms are defined by an iterated sum representation: 

x' 1 x' k 

Li mi> ,., j m k (xi,...,xj c ) = 52 7~mT "'Jm' ^ 

h>i 2 >...>i k >0 11 lJc 

The order of the arguments of multiple polylogarithms (and of multiple zeta values) is re- 
versed with respect to the definitions in [1,36-38]. It is also convenient to define the functions 
G(zi,—,Zk',y) for Zk 7^ by an integral representation as follows [3,38]: 




o o o 



Note that in the definition of the functions G(zi,...,Zk',y) one variable is redundant due to the 
following scaling relation: 

G(zi,---,Zk;y) = G(xzi,...,xzk\xy) (11) 

Within the functions G(zi,-—,Zk',y) one allows in addition trailing zeros (e.g Zk = Zk-l = ■■■ = 
Zj = 0) through the definitions 

G((W);y) = ^(lny)*, 

k 

y 

G(zi,...,Zk;y) = / G(z2,---,Zk\t). (12) 

J t-zi 
o 

For Zk 7^ this coincides with the original definition eq. (flOb . Note that the scaling relation (fTTT) 
does not hold for trailing zeroes. G-functions with trailing zeroes can always be transformed to 
G-functions without trailing zeroes. An algorithm for this transformation is given in section |4j 
Therefore we will always assume that Zk ^ 0, unless explicitly stated otherwise. The functions 
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G(zi,...,Zk',y) an d Limi,...,m k (xi, ■■■■>Xk) denote the same class of functions. With the short-hand 
notation 

G mu ..., mk (zi,...,zk;y) = G(0^^ 7 zu--^zk-\,0^0,zk-,y) (13) 

m\—\ mk~^ 

the relation between the two notations is given by 

Li mi; ..., mjt (xi, ...,Xk) = (— 1) Gmi,...,m t ( — i— — ■>■■■■>- ~T ; M" (l^) 

\X\ X\X2 X\...Xk J 

The notation Li mi) ... jinjt (;ti, ...,Xk) has a closer relation to the series representation, whereas the 
notation G(zi,---,Zk',y) has a closer relation to the integral representation. For the numerical eval- 
uation both representations will be important. It is further convenient to introduce the following 
notation for iterated integrals: 

A A t n -2 t„-\ 

dt dt f dt\ f dt n -i f dt n 

0...0 = / x ... x / / . (15) 

t — a\ t — a n J t\—a\ J t n -\—a n -\ J t n — a n 



We further use the following short hand notation: 

dt \ m dt f dt dt dt 

— o = / — o... — o . (16) 

t J t — a J ^£ ^ t y t — a 

m times 

With this notation, the integral representation of the function G m ^_^ mk (z\,Zi, ...,Zk',y) is written 
as 

}fdt\ mi - 1 dt (dt\ mi - X dt (dt\ mk - 1 dt 

G m ,...,m k {z\,Z2, ...,Zk\y) = )[jo) —^J*) T^ 2 -\T) — k 

(17) 

The series expansion for Li miv .. )mi (jci, is convergent, if 

|xiX2...xy| < 1 for all j e {1, ...,&} and (m\,x\) ^ (1, 1). (18) 

Therefore the function G m , mk (zi, ...,Zk',y) has a convergent series representation if 

\y\<\zj\ for all j, (19) 

e.g. no element in the set {|zi|, |zjt|, \y\} is smaller than \y\ and in addition if mi = 1 we have 
y/z\ 7^ 1. The power series expansion for the G-functions reads 

G mu ... 7 m k (z\,...,zk,y) = 

^i'"£i(7i + -+A) Wl (*) ( ./:•••••. /a r (|) -(AT© ' (20) 



5 



Multiple poly logarithms satisfy two Hopf algebras. The first one is referred to as "shuffle alge- 
bra" and is related to the integral representation. An example for the multiplication is: 



G(zi;y)G(zr,y) = G(zi,zr,y) + G(z2,zv,y). (21) 

The second algebra is often called "stuffle algebra" or "quasi- shuffle algebra" [39] and is related 
to the series representation. An example for this multiplication is: 

(JC1X2). (22) 

Special cases 

Multiple polylogarithms contain several specific subclasses. We list here the most important 
ones. The notation for multiple zeta values [2] is: 

fi>; 2 >->/*>0 11 ik 
Harmonic polylogarithms [5] are denoted as 

^Wl, (-*■) Lj-mj , ...,mf, ip^i 1 ; • • •1 1 ) ■ (24) 

k-1 

Nielsen's polylogarithms [7] are denoted as 

S n , P {x) = Li n+ i ) i r .. j i(x,l^T). (25) 

p-i 

Obviously, the class of multiple polylogarithms contains also the classical polylogarithms [6]: 

00 J 

3.2 Definitions related to the analytic continuation 

The multiple polylogarithms are analytic functions in k complex variables. To discuss the branch 
cuts it is convenient to consider the integral representation G mu ,„ : , nk (zi,Z2, ■■■,Zk\y) with Zk 7^ 0. 



tfdt r 1 " 1 dt fdt \ m2 - 1 dt fdt \ mk ~ ] dt 

G mi ,..., mk (zuz 2 ,..,z k ;y) = J (^oj — (^oj J^-\ T °) 



t-z k 
(27) 



Using the scaling relation eq. (fTTT) we can ensure that y is a positive real number. (In fact one 
could even require that y = 1, but y positive and real is sufficient for our purpose here.) The Zj are 
arbitrary complex variables. Therefore we can deduce from the integral representation, that the 
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function G mu ..^ mk (zi,Z2, ■■■,Zk',y) develops a branch cut whenever Zj is a positive real number and 
the integration variable t exceeds Zj. In most physical applications, the Zj will be real numbers. 
To distinguish if the integration contour runs above or below a cut, we define the abbreviations 
z±, meaning that a small positive, respectively negative imaginary part is to be added to the value 
of the variable: 

z+ = z + iO, z-=z-iO. (28) 

Close to the real axis we have 

Im l —— = ±%^-®(t-z). (29) 

t-zTiO dt 

This formula can be used to extract the imaginary parts from the G-functions [20]. 



4 Useful properties of multiple polylogarithms 

In this section we give an algorithm for the removal of trailing zeroes from G-functions, discuss 
the transformation properties of G-functions G(zi, ...,Zk',y) with respect to the argument y, inves- 
tigate in more detail the Bernoulli transformation and discuss the Holder convolution. The last 
two transformations can be used to speed up the series expansion. 



4.1 Trailing zeroes 

In this subsection we give an algorithm to remove trailing zeroes from G-functions. The method 
has been described for harmonic polylogarithms in [5] and has a straightforward generalisation 
to multiple polylogarithms. We say that a multiple poly logarithm of the form 

G(zi,...,zy,<W);;y) (30) 
k-j 

with zj 7^ has (k — j) trailing zeroes. Polylogarithms with trailing zeroes do not have a Taylor 
expansion in y, but develop a logarithmic singularity in (lny) around y = 0. In removing the 
trailing zeroes, one explicitly separates these logarithmic terms, such that the rest has a regular 
expansion around y = 0. The starting point is the shuffle relation 

G(0;y)G(z 1 ,...,z J -,0 i:::i 0;y)= (31) 

k-j-i 

(k-j)G(zh-,Zj,0,...,0;y)+ £ G(s h ...,Sj,Zj,0,...,0;y). 

{s h ...,sj)=(zi,...,Zj-i)UJ{0) ^ 

(zi, ...,Zj-i) LU (0) denotes the shuffle product of the string (zi, with (0). Solving this 

equation for G(zi,...,zy,0, ...,0;y) yields 

(32) 
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k-j 



G(0;y)G(zi,...,Zj,0,...,0;y)- £ G(s h ...,Sj,Zj,0,...,0;y) 

(s l ,...,Sj) = (z 1 ,-,Zj-i)U-l(0) 



In the first term, one logarithm has been explicitly factored out: 

G(0;y) = lny. (33) 

All remaining terms have at most (k — j — 1) trailing zeroes. Using recursion, we may therefore 
eliminate all trailing zeroes. 

4.2 Transformation of the arguments 

In the following we will consider transformations of 

G(z h ...,z k ;f(y)) (34) 

with respect to the argument y. For f(y) = 1 — y, f(y) = 1/(1 —y), f(y) = 1/y and f(y) = 
(1 — y)/(l +y) we give algorithms, which transform the G-function back to G(zi, ...,z^;y). For 
harmonic poly logarithms the corresponding transformations have been discussed in ref. [5]. The 
algorithms discussed here are straightforward extensions. 

The transformation 1 — y 

The function G(zi, ...,Zjt; 1 — y) can be expressed in terms of G-functions with argument y and 1 
as follows: 

y 

G{z\±,...,z k ;l-y) = G(zi ± ,...,^;l)+ / - — — G(z 2 , ...,z*; 1 • (35) 

o + 

For the second term we may use recursion. 
The transformation 1/(1 — y) 

The function G(zi, ...,Zk', 1/(1 —y)) can be expressed in terms of G-functions with argument y 
and 1 as follows: For zi ^Owe have 

G ^zi±,...,zfc; Y^^j =G(zi ± ,...,Zjfc;l) 



G ( Z2 ' rb) " / ^TI G ( Z2 ' T^) • 



(36) 
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For zi = we have 



G | 0,0, „.,0jZ2±, 1 = G | 0,0,...,0,Z2±,— ,Zjt;l 

mi— 1 / \ mi— 1 

>' 

dt I „ 1 



| ' G Lo, Z2± ,..., a; _ . (37) 



^ \ mi -2 

The transformation 1 jy 

The function G(zi, ...,Zk', 1/y) can be expressed in terms of G-functions with argument y and 1 
as follows: For z\ ^ we have 

1 y 

1\ ~, f dt „ / 1\ f dt „ / 1 



(38) 



G( zi±,...,Zjt;- J =G(zi±,...,Zjt;l) + y yG I z 2 , y J -y yG ( zi, ...,zjt; y 
37 o o 

i y 

/ * G ( Z2 ,..., a; I) + / * G ( Z2 ,..., a; I). 



For zi = we have 



g I o,o,...,o,z 2 ±,---,Z/t;- I =g I o,o,..., o, Z2±,— 



m\—\ / \ m\—\ 

•/>(—- /f-' 



0,...,0,Z2±,...,z^;y | • (39) 



mi -2 / \ mi-2 

Note that the 1/y transformation can also be obtained as a 1/(1 — y) transformation followed by 
a 1 — y transformation. 

The transformation (1 — y)/(l +y) 

The function G(zi , . . . , Zk', ( 1 — y) / ( 1 +y) ) can be expressed in terms of G-functions with argument 
y and 1 as follows: For zi ^ — 1 we have 

G^zi±,-",Z/byy^ = G(zi ± , -,Zjt; 1) 

+ / ck,...,*!^- /^G^,...,^. (40) 



f / 1-Zl \ \ 1+tJ J t+l \ l+t 

1 \\+zx J T 
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For zi = — 1 we have 




4.3 Bernoulli substitution 

In this subsection we investigate in more detail the Bernoulli transformation, which can be used 
to speed up the series expansion. We recall the standard series representation for the dilogarithm: 

Li 2 (x) = (42) 
j=i J 

Due to eq, © and eq. © we can assume that |x| < 1 and Re(x) < 1/2. Eq. (|42b converges 
therefore geometrically with x. On the other hand: 



x -ln(l-x) 

Li2W = -f± ln(l- t )= / ^=£.77^ 





where we used 



z J+1 , (43) 



it « — 1 t/l 

= YBj— (44) 

« 1 7=0 

and the notation z = — ln(l — x). For large j one obtains from the formula 

bj = (- 1 ) (i+2)/2 (|^7^' ;' even > ( 45 ) 

that the sum converges geometrically with zj (2%), e.g. 

5 



(7 + 1)! 

For x in the range [—1,1/2] we have 



In(l-x) 



2% 



l)^V*fl) J . (46) 



< — < 0.1104. (47) 
2n 



Therefore we can expect an improvement in the speed of convergence by using eq. d4*3l instead 
of eq. (l4*2"l) . We can generalise the Bernoulli transformation to the classical polylogarithms and 
to the harmonic polylogarithms. For the classical polylogarithms we find 

Li «« = l7^1lTT(- ln (l-^)) J+1 , (48) 
10 



where C\ (j) = 8 J; o and 

c - +iW = £(0^ c " w - <49) 

The coefficients C n (j) are independent of x and need therefore to be calculated only once. Similar 
we obtain for the harmonic polylogarithms: 

H m ,..,m k (x) = £^f^(-ln(l-^)y +1 , (50) 

where 

Cl ' m2 ' - ^ a) = { Cm 2 ,...,m k (j-l), J j>0\ (51) 

and 

C mi + l,m 2 ,...,m k U) = 52 ^ k ^ j"Cni,m 2 ,...,m t (^)- (52) 

For the classical polylogarithms hi n (x) the Bernoulli transformation leads to a speed-up if n is 
not too large. Empirically we find that for n > 12 the "standard" series expansion without the 
Bernoulli transformation is faster. This can be related to the explicit factor 1 / j n in the series 
expansion 

oo j 



4.4 Holder convolution 



The multiple polylogarithms satisfy the Holder convolution [2]. For zi ^ 1 and z w 7^ this 
identity reads 

G(zi,...,z w ;l) = J^(-l)- 7 G^l-z;,l-z;_i,...,l-zi;l-^ G (zj+i, ...,z w ; ^ . 



(54) 

The Holder convolution can be used to improve the rate of convergence for the series represen- 
tation of multiple polylogarithms. 
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5 Numerical evaluations 



5.1 Numerical evaluation of multiple zeta values 

Every multiple zeta value C,(s\, . ..,Sk) with s\ > 1 is finite and can be written as a convergent 
sum. Therefore it seems that the numerical evaluation is straightforward and no further trans- 
formations are needed. Yet the convergence can be quite slow if the parameters are small or 
high precision is needed. Therefore algorithms have been invented to accelerate the evalua- 
tion [31,40,41]. The motivation often was to find numerically new relations among multiple zeta 
values, and to this end the algorithms needed to deliver a very high digit accuracy, i.e. several 
hundred decimal digits, efficiently in short time. Of course, low precision evaluations benefit 
from these algorithms as well and since multiple zeta values with small parameters appear in the 
transformations of harmonic and multiple poly logarithms, these acceleration techniques are also 
important for the numerical evaluation of those functions. 

The algorithm proposed by Crandall [31] partitions the iterated integral by a chosen constant 
X, in such a way that the multiple zeta values can be expressed as a finite sum over two iterated 
integral functions Y and Z, for which all integration variables are either smaller or greater than X 
respectively. Different acceleration method apply for Y and Z. 

Z can be written in almost exactly the way as C, but has an additional factor e~^ x {n\ is the 
outermost summation index). This factor improves the convergence, and the greater X is the 
faster Z converges. 

If X is chosen to be smaller than 2n, Y can be rewritten in terms of Bernoulli numbers. The 
remaining summation can be interpreted as a convolution and this convolution can be made fast 
by applying FFT methods. The vectors to be convoluted and containing Bernoulli numbers or 
factors of gamma functions have to have a certain length L, which becomes a parameter of the 
algorithm like %. The performance of the evaluation depends on the values of X and L, but quite 
generally one can find that Crandall's algorithm performs best for high precision and a small 
number of C, parameters. 

Another method [2] uses the Holder convolution. As a generalisation of duality a given 
£>mi,...,m k is rewritten as a finite sum over functions 

G Wl) ... jmjk (l,...,l;l/p) and G mi) ... jmjfc (1,...,1;1/?). (55) 

For p,q>\ they have a significantly better convergence than £ mi; ... ; m r The p and p must satisfy 
the Holder condition l/p + l/q = 1 and can be chosen for best convergence as p = q = 2. 

The acceleration is significant and exceeds that of Crandall's algorithm for low precision 
evaluations. Both algorithms can play hand in hand to satisfy different demands. For high preci- 
sion evaluations Crandall's algorithm should be chosen, for the other cases Holder convolution 
is the best choice. 

5.2 Numerical evaluation of harmonic polylogarithms 

Harmonic polylogarithms have been introduced by Remiddi and Vermaseren [5] and their numer- 
ical evaluation has been worked out by Gehrmann and Remiddi [32]. Harmonic polylogarithms 
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have a convergent series representation if \x\ < 1 and no trailing zeroes are present. The right- 
most indices are allowed to be zero, in which case the harmonic polylogarithm does not have a 
series expansion around zero in x, but develops a logarithmic singularity proportional to lnx To 
cope with this trailing zeros scenario, one uses recursively the shuffle algebra to extract those 
zeros in form of ordinary logarithms. It remains to transform an arbitrary argument x into the 
domain \x\ < 1. To this aim one uses the transformations I fx and (1 — x)/( \ +x) from sect. 
14.21 Here we differ slightly from the implementation by Gehrmann and Remiddi: These authors 
provide a numerical evaluation routine for harmonic polylogarithms up to weight 4. They use 
for x G [Vz — 1,1] the ( 1 — x)/{ 1 + x) transformation to improve the convergence of the result- 
ing series. We provide a numerical evaluation routine for arbitrary weights and therefore have 
to determine the transformation formula at run time. We observed that this transformation by 
itself is quite expensive. This can outweigh the potential gain from the series acceleration. It is 
therefore more economical to apply this transformation only for x near one. The starting point is 
determined empirically. 

In addition, multiple zeta values appear in the transformed expressions. If some indices 
are negative, alternating multiple zeta values occur. The speed of the numerical evaluation of 
multiple zeta values also has a major impact on evaluation speed of the harmonic polylogarithms. 



5.3 Numerical evaluation of multiple polylogarithms 

Multiple polylogarithms can be evaluated with the help of the series expansion in the domain 
where this expansion is convergent. This is the case if the arguments satisfy eq. (fT8l) or equiv- 
alently eq. dT^l. We first show that any combination of arguments can be transformed into 
this domain. This transformation is based on the integral representation and is therefore most 
naturally expressed in terms of the G-f unctions. We consider the function 



Gmi,...,mj (zi,...,Zj-l,S,Zj+l, ...,Zk',y) 



(56) 



with the condition that \s\ is the smallest element in the set {|zi|, \s\, \zj+i\, \zk\, \y\}- 
The algorithm goes by induction and introduces the more general structure 



ds] 



si — b\ 



dSr 



s r -b r 



G(a\,...,s r , ...,a w ;y2), 



(57) 



where |yi| is the smallest element in the set {|yi|, \b\\, \b' r \, \a[\, \a' w \, \yi\}- The prime in- 
dicates that only the non-zero elements of at and bj are considered. If the integrals over si to s r 
are absent, we recover the original G-function in eq. d56l) . Since we can always remove trailing 
zeroes with the help of the algorithm in section l4~Tl we can assume that a w ^ 0. We first consider 
the case where the G-function is of depth one, e.g. 



ds] 



S r -1 



dSr 



y\ 



-G(p,...,0,s r ;y 2 ) 



ds\ 



s r -\ 



m—1 



si -b\ 



dSr 



-G m (s r ;y 2 ), (58) 
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and show that we can relate the function G m (s r ;y2) to G m (y2',s r ), powers of ln(sy) and functions, 
which do not depend on s r . For m = 1 we have 

Gi(s r ±;n) = G 1 (y 2T ;s r )-G(0;s r )+ln(-y2 T ). (59) 
For m>2 one can use the transformation 1 jy from sect. I4.2l and one obtains: 



G m (s r± ;y 2 ) = -^ m + J -jG m -i(t±\y2)- J -^G m -\(t±\y2). 



(60) 



One sees that the first and second term in eq. (l60l) yield functions independent of s r . The third 
term has a reduced weight and we may therefore use recursion. This completes the discussion 
for G m (s r ;y2). We now turn to the general case with a G-f unction of depth greater than one in 
eq. (l57t . Here we first consider the sub-case, that s r appears in the last place in the parameter list 
and (m — 1) zeroes precede s r , e.g. 

/ r~ / r -rG(ai,...,a k ,0,...,0,s r ;y2), (61) 

J S]-b] J s r -b r ^ — „ — ' 



m—1 



Since we assumed that the G-function has a depth greater than one, we may conclude that ^ 0. 
Here we use the shuffle relation to relate this case to the case where s r does not appear in the last 
place: 

G(ai,...,%,0,...,0,^;y2) = G(«i, ...,a k ;y 2 )G(0, ...,0,s r ;y 2 ) - G(ai,...,a yt+m ;v 2 ), 

V "T v i shuffles' 

m— 1 m— 1 J J 

(62) 

where the sum runs over all shuffles of (ai, ...,afc) with (0, ...,0,s r ) and the prime indicates that 
(oci, ...,a^ +m ) = (ai,...,%,0, ...,0,* r ) is to be excluded from this sum. In the first term on the 
r.h.s of eq. d62l the factor G(ai, 72) is independent of s r , whereas the second factor 
G(0, 0, s r ; V2) is of depth one and can be treated with the methods discussed above. The terms 
corresponding to the sum over the shuffles in eq. (l6"2"l) have either s r not appearing in the last 
place in the parameter list or a reduced number of zeroes preceding s r . In the last case we may 
use recursion to remove s r from the last place in the parameter list. It remains to discuss the case, 
where the G-function has depth greater than one and s r does not appear in the last place in the 
parameter list, e.g. 

y\ Sr-i 

/ \- / r —G(ai,...,ai-i,s r ,ai + i,...,a w ;y2), (63) 

J s\—b\ J Sy — Or 


with a w 7^ 0. Obviously, we have 

G(ai,...,ai-i,s r ,ai + i,...,a w ;y 2 ) = (64) 

? d 

G(ai, ...,ai-\,0,ai + \,...,a w ;y2) + / ds r+ \- G(ai, ...,a,-_i, j r+ i,a,-+i, ...,a w ;y 2 ) ■ 

J ds r+ i 



14 



The first term G(a\, ...,a;_i,0,a/+i, ...,a w ;y 2 ) does no longer depend on s r and has a reduced 
depth. For the second term we first write out the integral representation of the G-function. We 
then use 

ost — s at t — s 

followed by partial integration in t and finally partial fraction decomposition according to 

11 1/1 1 



t — at — s s-a \t — s t-a 
If s r is not in the first place of the parameter list, we obtain 



(66) 



[ ds r+ \zr G(ai,...,ai-i,s r +i,ai+i,...,a w ;y2) 

J ds r+ i 
o 

Sr 

= — / — — G(ai, ...,di-2,s r+ i,ai + i, ...,a w ;y 2 ) 



o 



+ / G(ai,...,aj_2,ai_i,ai + i, ...,a w ;y2) 

J S r+ l—ai-l 


+ / G(ai, ...,ai-i,s r+ i,ai+2,...,a w ;y2) 

J s r+ \-a i+ i 
o 

Sr 

- / — — G(ai,...,a l -_i,a,- + i,a l - + 2,...,a w ;y2)- (67) 







Each G-function has a weight reduced by one unit and we may use recursion. If s r appears in the 
first place we have the following special case: 



Sr _ s r 



[ ds r+ i- G(s r +i,ai+\,...,a w ;y2) = / G(a ;+ i, ...,a w ;j2) 



+ / — — G(s r+ i,ai + 2,...,a w ;y 2 )- f — — G(a i+u a i+2 , ...,a w ;y 2 ) ■ (68) 
J s r+ \—ai + \ J s r+ \-a i+ \ 

o o 

There is however a subtlety: If a,-_i or a !+ i are zero, the algorithm generates terms of the form 

j^F(s)-j^F(0). (69) 



o o 
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Although the sum of these two terms is finite, individual pieces diverge at s = 0. We regularise 
the individual contributions with a lower cut-off X: 




In individual contributions we therefore obtain at the end of the day powers of InX from integrals 
of the form 

j d JlJ d Jl = Wy-lnylnl + Wl. (71) 
J s\ J S2 2 2 

X X 

In the final result, all powers of InX cancel, and we are left with G-functions with trailing ze- 
ros. These are then converted by standard algorithms to G-functions without trailing zeros. The 
G-functions without trailing zeros can then be evaluated numerically by their power series ex- 
pansion. 

In addition, the algorithms may introduce in intermediate steps G-functions with leading 
ones, e.g G(l, ...,z&;l). These functions are divergent, but the divergence can be factorised and 
expressed in terms of the basic divergence G(l; 1). The algorithm is very similar to the one for 
the extraction of trailing zeroes. In the end all divergences cancel. 



Acceleration of the convergent series 

The G-function G mi) ... jmjt (zi, —,Zk',y) has a convergent power series expansion if the conditions 
in eq. (IT9t are met. This does not necessarily imply, that the convergence is sufficiently fast, 
such that the power series expansion can be used in a straightforward way. In particular, if z\ 
is close to y the convergence is rather poor. In this paragraph we consider methods to improve 
the convergence. We can assume that no trailing zeroes are present (z% ^ 0), therefore we can 
normalise y to one. Convergence implies then, that we have \zj\ > 1 and (zi,mi) ^ (1, 1). If 
some Zj is close to the unit circle, say, 

1<M<2, (72) 
we use the Holder convolution eq. d5?b with p = 2 to rewrite the G-functions as 

G(zi,...,z w ;l) = G(2 Z i,...,2 Zw ;l) + (-l) w G(2(l-z w ),2(l-z w _i),...,2(l- Zl );l) 

w— 1 

+ £ (-1)^(2(1-2,0, 2(1 -Zj-i),...,2(l-zi);l)G(2z j+ i,...,2z w ;l). (73) 
;=i 

Here, we normalised the r.h.s to one and explicitly wrote the first and last term of the sum. 
We observe, that the first term G(2zi, ...,2z w ; 1) has all arguments outside \2zj\ > 2. This term 
has therefore a better convergence. Let us now turn to the second term in eq, d73l . If some Zj 
lies within \zj — 1| < 1/2, the Holder convolution transforms the arguments out of the region 
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of convergence. In this case, we repeat the steps above, e.g. transformation into the region of 
convergence, followed by a Holder convolution, if necessary. While this is a rather simple recipe 
to implement into a computer program, it is rather tricky to proof that this procedure does not lead 
to an infinite recursion, and besides that, does indeed lead to an improvement in the convergence. 
For the proof we have to understand how the algorithms for the transformation into the region 
of convergence act on the arguments of a G-function with length w. In particular we have to 
understand how in the result the G-functions of length w are related to the original G-function. 
Products of G-functions of lower length are "simpler" and not relevant for the argument here. We 
observe, that this algorithm for the G-function G(zi,...,z w ;y) substitutes y by the element with 
the smallest non-zero modulus from the set {\zi | , \z w \ , \y\}, permutes the remaining elements 
into an order, which is of no relevance here and possibly substitutes some non-zero elements by 
zero. The essential point is, that it does not introduce any non-trivial new arguments (e.g. new 
non-zero arguments). 

For the Holder convolution we are concerned with the second term of eq. (1731) 

G(2(l -z w ),2(l - Zw -i),...,2(l-zi); 1). (74) 

The first term G(2zi, ...,2z w ; 1) never transforms the arguments into the non-convergent region 
and has, as we have seen, a better convergence. The terms in the sum of eq. (1731) have a reduced 
length, and by induction we can assume that suitable methods to improve the convergence exist 
for those terms. For the second term of eq. ( f73t we have to discuss the transformation z — > 
2(1 — z). We divide the arguments Zj into different classes: 

- class A: \z\ > 2. These will map under the transformation z — > 2(1 — z) again to \z\ > 2. 

Actually, they transform to \z — 2| > 4, but this region is included in the previous one. 

- class B: 1 < |z| < 2 and \z — 1 1 > 1 . These will map under the transformation z — > 2( 1 — z) to 

\z\ > 2 (and necessarily \z — 1| > 1). Therefore class B is mapped to class A. 

- class C: 1 < |z| < 2 and 1/2 < |z — 1| < 1. These will map under the transformation z — > 

2(1 — z) to 1 < |z| < 2 and \z — 1| > 1. Therefore class C is mapped to class B. 

- class D: 1 < |z| < 2 and \z — 1| < 1/2. These will map under the transformation z — > 2(1 — z) 

to \z\ < 1 and \z — 1| > 1. Class D is mapped to the non-convergent region \z\ < 1. 

Let us first assume that all z/s are from class A and B. Then after the Holder convolution, all 
arguments satisfy |z| > 2, which ensures a fast convergence. Secondly, we assume that all z/s 
are from the classes A, B and C. The the Holder convolution will generate G-functions with 
arguments from the classes A and B alone. One subsequent Holder convolution on those G- 
functions, which contain arguments from class B will again lead to \z\ > 2 for all Zj- For the 
last case we have to consider G-functions, where some arguments are from class D. Then we 
obtain arguments with |z| < 1 and it is necessary to re-use the transformation into the region of 
convergence. Let nco be the number of arguments Zj from the classes C and D from the original 
G-function. The second term in eq. (1731) equals eq. d74"l) . Let z m in be the argument such that 
2 1 1 —z m in\ is the smallest in the set 

{2|l-z w |,..,2|l-z 1 |,l}. (75) 
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Since at least one argument is from class D, we have 2\1 — z m in\ < 1- The algorithm for the 
transformation into the convergent region introduces then G-functions of the form 



1 



7+1 



1 



1 -Zo,-_, 



1 - Zoi . 



(76) 



1 Zmin 1 Zmin 2(1 Zmin) 1 Zmin 1 -Zmin 

Here, a is permutation of the original arguments. The argument 1/2/(1 — Zmin) results from 
permutating the original y = 1 into the argument list of the z's. We note that this argument 
satisfies 

1 



2(1 



1 



Zmin J 



>1- 



(77) 



and lies therefore always outside region C and D. This is most easily seen by using polar coor- 
dinates for 2(1 — Zmin)- Furthermore, if an arbitrary argument z is from class A or B, it satisfies 
\z— 1| > 1. Therefore 



> 



1 



and it follows that 



1 



1 



" Zmin 



> 1. 



(78) 



(79) 



Therefore, arguments from classes A and B will remain in these classes. It remains to consider 
the case z = and to show it does not re-introduce arguments in the classes C or D. Again, we 
can show that 



1 



1 



1 



> 1. 



(80) 



In summary, the G-function in eq. (I76T) has the number hcd reduced by at least one (due to eq. 
1771 . Therefore the algorithm will successively remove the arguments from classes C and D and 
terminate. This completes the proof. 

In practice, there is again a trade-off between the gain in acceleration and the cost involved 
for the Holder convolution. We therefore apply the Holder convolution only if some zj satisfies 

l<\zj\<X, (81) 



with X < 2. A typical value is X = 0.01. 



6 Implementation and checks 

The numerical evaluations have been implemented as part of GiNaC [35], a C++ library for 
computer algebra Qhttp : / /www . qinac . de| ). GiNaC enables symbolic algebraic manipulations 
within the C++ programming language. Like FORM [42], it was developed within the high- 
energy physics community. GiNaC allows numerics in arbitrary precision. 

The following functions have been added to GiNaC. On the left are the names and the math- 
ematical notations, on the right are the names within GiNaC. 
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Name 



Math, notation GiNaC notation 



Classical polylogarithm 
Nielsen polylogarithm 
Harmonic polylogarithm 
Multiple polylogarithm 



Li„(x) Li(n,x) 

S n:P (x) S(n,p,x) 

Ufft(x) H({ml, . . .,mk},x) 

Li^(x) Li ( {ml, . . . ,mk} , {xl, . . . , xk} ) 

G(a,y) G({al, . . .,ak},y) 



G({al, . . .,ak}, {dl, . . .,dk},y) 



Multiple zeta value 
Alternating MZV 



C,(m) zeta ( {ml, . . . , mk} ) 

^(m,d) zeta ( {ml, . . . , mk} , {si, . . . , sk} ) 



While x and the entries of x can be arbitrary complex numbers, n, p and the entries of m 
must be positive integers. In the case of the harmonic polylogarithm m may also contain negative 
integers, a is the expanded parameter string for the G-functions. The dfs are positive or negative 
numbers indicating the signs of a small imaginary part of a,-. The s ; 's are positive or negative 
numbers indicating which corresponding sum in the definition of the multiple zeta values is 
alternating. 

The algorithms used can be found in [28, 29] for classical polylogarithms and Nielsen poly- 
logarithms, in [5, 32] for harmonic polylogarithms, in [2, 31] for multiple zeta values and in 
our text above for multiple polylogarithms. The transformations are performed algebraically so 
there are no restrictions on the size or length of the parameters. The summation is accelerated by 
the use of Bernoulli numbers in the case of classical and Nielsen polylogarithms. Since GiNaC 
has arbitrary precision numerics and the precision can be changed at every time this is the only 
feasible acceleration technique. Evaluation of multiple zeta values either goes by Crandall's al- 
gorithm or by Holder convolution depending on the precision and parameters. The switching 
point has been determined empirically. For alternating multiple zeta values Holder convolution 
is always used. 

For the sign of the imaginary part we adopted a widely used convention for mathematical 
software: implementations shall map a cut so the function is continuous as the cut is approached 
coming around the finite endpoint of the cut in a counter clockwise direction [43]. With this 
convention the cuts on the positive real axis are continuous to the lower complex half-plane. 
This convention also ensures consistency among the various polylogarithms including the ordi- 
nary logarithm. We note that the implementation of harmonic polylogarithms by Gehrmann and 
Remiddi [32] uses the opposite sign convention. Gehrmann and Remiddi give the argument of H 
a positive imaginary part x + ie. This implies that cuts on the positive real axis are continuous to 
the upper complex half -plane. 

To ensure the correctness of the implementation many checks have been performed. First, 
one can compare the numerical evaluation to known special values of the function, e.g. Li„ 
with parameters x = 0, \ or 1, or ^(2n) with integer n. In the same sense but with less rigour of 
course, one can compare to other well proven software implementations for a quick reassurance. 
Second, the transformations can be checked for internal consistency if they overlap. Hence the 
\—x transformation can always be be verified against the pure series. If the (1 — x)/{\ +x) 
transformation is done as it is with the harmonic polylogarithm the l/x transformation overlaps 
and can be tested as well. Checking identities between different functions is the third method. 
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Since all discussed functions merely are special cases of the general multiple polylogarithm, 
there exists a hierarchy of identities between them, e.g. 

S n ,i(x) = Li M+ i(x) or H 3i i 5 i(x) = S 2 ,3(»- 

In addition, limits can be examined, for example 

Li2,3,4 (x,a, 1) -> H 2i 3 i4 (x) as a -> 1. 

The checks listed above are sufficient to assure the correctness of the implementation for the 
simpler functions. Yet the more general functions of multiple poly logarithms and multiple zeta 
values still have areas of the parameter space that have not been tested with the methods so far. 
Here, the shuffle and quasi-shuffle identities are a very important tool to fill this gap. Identities 
like 

C(4,3) = 17C(7)-10C(2)C(5) 

or 

Li 2 (x)Li 5 (y) = Li 2j5 (x,y) + Li 5;2 (;y,x) +Li 7 (xy) 

are examples for such identities. Shuffle and quasi-shuffle relations similar to those can be used 
to expose problems within the implementation for arbitrary parameter sets. 

The checks mentioned above have been exercised not only once, but most of them have been 
added as automatic checks to the GiNaC library itself. With the source code at hand the test suite 
of GiNaC can be run anytime. This helps to detect new errors quickly that might make their way 
into the code in the course of future extensions or optimisations. 

The standard use of GiNaC is through its C++ interface. A small example program is listed 
below: 

#include <iostream> 
#include <ginac/ginac . h> 

int main() 

{ 

using namespace std; 
using namespace GiNaC; 

ex xl = numeric (8, 3) ; 
ex x2 = numeric (1, 5) ; 

cout << "Li_{l, 1} (8/3, 1/5) = " 

<< Li( lst(l,l), lst(xl,x2) ).evalf() « endl; 

return 0; 

} 

If this program is compiled and linked with the GiNaC-library, it will print out: 
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Li_{ 1,1} (8/3,1/5) = -0.8205920210842043836-0.70102614150465842094*1 

Alternatively, GiNaC offers also a small interactive shell called ginsh, which allows to try and 
use GiNaC's features directly as in the following examples: 

> Li (2,1); 
l/6*Pi A 2 

> S (2,3,4.5) ; 

-1 .52140580215075747 68 + 1 .7 01377 689228 92 6854 6*1 

> Li({2, 2,1}, {3.0,2.0,0.2}); 

-0.7 89067882 6631402472+0.57 916837032172 81085*1 

> Digits=40; 
40 

> H({2,-1,3},8.7); 

-5.652 07410697321998445159060623787475178342 96803 6-1 .0548 62 93307539105482 
502537 832457314244 07 02785858*1 

Here the user input is done at the prompt > and the result is given on the next line. 



7 Conclusions 

In this paper we reported on numerical evaluation methods for multiple polylogarithms. These 
functions occur in higher loop calculations in quantum field theory. We provided algorithms, 
which allow the evaluation of these functions for arbitrary complex arguments without any re- 
striction on the weight of the function. The functions can be evaluated in C++ to arbitrary 
precision within the GiNaC framework. Subclasses of multiple polylogarithms are the classical 
polylogarithms, the Nielsen polylogarithms, the harmonic polylogarithms and the multiple zeta 
values. For these subclasses we have implemented specialised algorithms. All routines are in- 
tegrated in the computer algebra package GiNaC from version 1.3 onwards and can be obtained 
by downloading this library [35]. 
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